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HIGHLIGHTS 


•  1-d  axisymmetric  plane  strain  model  to  study  stress  in  different  Si  nanostructures. 

•  1-d  model  is  a  dimensionally  reduced  form  of  the  2-d  model  presented  in  part  I. 

•  1-d  and  2-d  model  results  compare  well  for  the  Si  nanowire. 

•  Si  nanotubes  expected  to  perform  better  in  terms  of  mechanical  stability. 

•  Carbon-Si  core/shell  structures  exhibit  high  stress  generation  at  interface. 
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The  second  part  of  the  manuscript  presents  a  one-dimensional  axisymmetric  plane  strain  model  to  study 
lithium  insertion  induced  stresses  in  different  types  of  silicon  nanostructures,  —  nanowire  (Si  NW), 
nanotube  (Si  NT)  and  core/shell  (Si  C/S)  nanostructures.  The  1-d  plane  strain  model  is  the  dimensionally 
reduced  form  of  the  2-d  model  presented  in  part  I.  Simulation  results  from  the  1-d  plane  strain  model  is 
compared  to  the  results  from  the  2-d  model  for  the  Si  NW.  Each  of  these  structures  poses  a  different  type 
of  boundary  conditions  and  the  stress  evolution  due  to  diffusion  in  each  of  the  cases  are  discussed. 
Simulation  results  indicate  that  Si  NT  structures  develop  lower  tensile  stresses  compared  to  Si  NW 
structures  under  similar  current  densities.  Case  studies  on  Si  NT  for  different  values  of  inner  to  outer  radii 
ratio  are  also  presented.  Furthermore,  simulations  reveal  that  the  Si  C/S  structures  develop  much  higher 
stresses  (closer  to  the  core/shell  interface)  compared  to  Si  NW  or  Si  NT,  owing  to  the  different  expansion 
factors  of  the  core  and  shell  material. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Nanostructured  electrodes  have  been  proposed  to  alleviate  the 
detrimental  effects  of  large  volume  changes  due  to  lithium  inser¬ 
tion  in  Si  electrodes.  Experiments  have  demonstrated  that  nano¬ 
structured  materials  [1-3]  can  better  accommodate  large  diffusion 
induced  stresses  caused  due  to  insertion  and  de-insertion  of 
lithium,  thereby  increasing  cyclability  of  Si  electrodes.  Many 
research  groups  have  reported  the  use  of  various  nano-structured 
materials  including  Si  nanowires,  Si  nanotubes,  core— shell  Si 
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nanowires  etc.  as  electrodes  in  lithium-ion  batteries  in  an  attempt 
to  increase  cycle  life  and  performance  of  the  materials  1-6].  Yao 
et  al.  [7  reported  interconnected  Si  hollow  nanospheres  as  Li-ion 
battery  anodes  with  high  cycle  life  as  a  result  of  reduced  stresses 
due  to  the  hollow  space  compared  to  solid  structures.  Song  et  al.  [4] 
shows  similar  reduction  of  stress  in  Si  NT.  Liu  et  al.  [8  reported  a 
yolk-shell  structure  consisting  of  Si  nanoparticle  enclosed  in  a 
carbon  shell  with  void  space  between  particle  and  shell  so  that 
particle  can  expand  freely  during  lithiation  without  breaking  the 
shell.  Core  shell  (carbon/Si)  nanostructures  have  also  been  reported 
in  literature  for  use  as  high  capacity  anodes  for  Li  ion  batteries  [3,9]. 
The  authors  claim  that  the  carbon  core  provides  mechanical  sup¬ 
port  to  the  nanostructure  system  making  it  more  stable  when  huge 
stresses  build  up  during  lithiation.  Other  novel  structures  have 
been  reported  in  literature  to  better  accommodate  volume  expan¬ 
sion  in  silicon,  such  as  a  carbon-aluminum-silicon  nanostructure 
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comprising  of  a  carbon  nanorod,  a  sandwiched  layer  of  aluminum, 
capped  by  a  nanoscoop  of  silicon,  [10  which  exploits  the  graded 
strain  properties  of  the  three  different  materials  to  accommodate 
the  volume  expansion  during  lithiation.  In  an  earlier  work  (Part  I  of 
this  paper),  we  presented  a  2-d  model  for  electrochemical  lithiation 
of  a  Si  nanowire.  In  this  article,  we  present  a  simplified,  1- 
dimensional,  computationally  efficient  mathematical  model  to 
simulate  stress  generation,  volume  expansion  etc.  due  to  lithium 
insertion.  The  1-d  model  is  used  to  simulate  lithiation  of  the 
following  silicon  based  nano-structured  electrodes,  -  nanowire  (Si 
NW),  nanotube  (Si  NT)  and  core/shell  (Si  C/S)  nanostructures. 

2.  Model  development 

A  detailed  transient  axisymmetric  two  dimensional  model  (2- 
d  model)  to  simulate  the  electrochemical  lithium  insertion  in  Si 
nanowire  has  been  developed  in  part  I  of  the  paper.  The  2-d  model 
was  reduced  to  a  1-d  model  in  r  co-ordinate  with  the  intent  of 
minimizing  computational  run  time  which  facilitates  faster  analysis 
of  multiple  structures.  The  key  assumptions  which  are  similar  for 
both  1-d  and  2-d  models  were  discussed  in  part  I.  Additional 
assumptions  involved  in  the  1-d  model  are  as  follows: 

•  The  total  strain  in  the  axial  z  direction  ( ezz )  is  assumed  to  be  zero 
since  the  axial  dimension  of  each  nanostructure  is  much  greater 
compared  to  the  radius.  This  assumption  neglects  the  variation 
of  all  variables  including  the  displacement  due  to  lithium 
insertion  in  the  z  direction  facilitating  the  development  of  a  one 
dimensional  axisymmetric  plane  strain  model.  The  1-d  plane 
strain  model  is  expected  to  be  representative  in  regions  away 
from  the  top  and  bottom  edges  of  the  nanostructures. 

•  Transport  of  Li  in  Si  host  due  to  diffusion,  pressure  gradient  and 
bulk  transport  is  considered  only  in  the  radial  direction  and 
transport  in  the  axial  and  tangential  co-ordinates  is  ignored. 

•  The  current  collector  constraint  at  the  base  of  nanostructure  is 
ignored  which  is  also  representative  of  region  far  away  from  the 
top  and  bottom  of  the  nanostructures. 

•  For  the  Si  NT,  the  inner  surface  is  assumed  to  be  devoid  of 
electrolyte,  thereby  ignoring  any  electrochemical  reaction. 

•  For  the  Si  shell  carbon  core  structures,  lithium  insertion  in  the 
carbon  core  is  neglected,  due  to  relatively  smaller  volume 
element  and  lower  expansion  of  carbon  compared  to  Si. 


shell  carbon  core  (Si  C/S)  nanostructures.  Fig.  1  shows  the  schematic 
of  the  structures.  The  1-d  model  is  treated  with  necessary  modifi¬ 
cations  for  each  nanostructure  along  with  appropriate  boundary 
conditions  to  simulate  each  case.  The  results  from  1-d  Si-NW  model 
are  compared  to  those  from  the  2-d  model. 


2.1.  Model  equations 


2AA.  Material  balance 

The  material  balance  equations  related  to  the  1-d  model  are 
similar  to  those  described  in  part  I  of  paper.  An  exception  is  that  for 
the  1-d  model,  axial  variation  is  neglected  and  quantities  vary  only 
across  r.  Therefore  using  similar  notations  as  in  part  I  of  the  paper, 
the  flux  of  lithiated  host  or  species  LiS  is  expressed  as  follows: 
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The  material  balance  for  species  LiS  is  given  by 
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2.1.2.  Solid  mechanics 

Due  to  the  plane  strain  assumption,  the  axial  strain  is  taken  to  be 
zero  i.e.  szz  =  0.  Therefore,  the  large  deformation  strain  matrix  is 
reduced  to  the  following  equation  for  the  1-d  case: 
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Here  u  is  the  radial  displacement  corresponding  to  elastic  strain. 
The  radial  and  tangential  stresses  are  related  to  the  strain  using 
linear  Hooke’s  Law  given  by 


t  =  Atr(<?)(I)+2  /i(e) 


(4) 


The  paper  presents  a  study  and  analysis  of  stress  distributions 
for  different  silicon  based  nanostructured  electrodes  using  the 
simplified  1-d  model.  The  different  nanostructures  considered  are 
a)  Silicon  nanowire  (Si  NW)  b)  Silicon  nanotube  (Si  NT)  c)  Silicon 


Case  1  Si  NW  2  D  Axisymmetric  Representation  1  D  Axisymmetric  Representation 
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Fig.  1.  Schematic  and  modeling  geometry  representation  of  Si  based  nanostructures 
considered  for  model  development. 


where  t  is  the  stress  tensor,  A  and  p  being  the  Lame  constants. 
The  corresponding  equilibrium  force  balance  equation  reduces 
to 

torr  ffrr-fftt  =  Q  (5) 

Or  r 

The  respective  one  dimensional  force  and  material  balances 
were  solved  in  a  similar  procedure  discussed  in  part  I  for  the  Si  NW, 
Si  NT  and  Si  shell  (for  core  shell  structures)  domains.  For  the  Si  shell 
carbon  core  structures,  the  model  variables  were  also  solved  inside 
the  carbon  core  domain.  The  equilibrium  force  balance  equation 
(Eq.  (5))  is  valid  in  the  core  domain.  As  per  our  assumption,  there  is 
no  Li  insertion  into  the  core,  so  the  total  strain,  eT  is  purely  elastic 
and  the  chemical  strain  component  is  absent  in  the  core.  Therefore 
the  radial  displacement  component  u  is  the  only  variable  to  be 
solved  for  in  the  core  domain. 

2.2.  Boundary  and  initial  conditions  for  different  nanostructures 
2.2 A.  Case  A.  Silicon  nanowire 

The  1-d  model  is  solved  for  a  constant  current  case.  It  should  be 
noted  that  the  outer  radius  specified  by  R0  changes  with  time  as  the 
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structure  expands  during  lithiation.  At  the  outer  surface,  the  cur¬ 
rent  flux  boundary  condition  for  constant  current  condition  is 
applied. 


-n-Nus 


\r= 0 


-Sjiapp(t) 

F 


(6) 


iapp  represents  the  current  density  which  is  a  function  of 
time  due  to  the  changing  outer  surface  area.  The  applied  current 
is  constant  which  satisfies  the  galvanostatic  conditions.  The 
axial  symmetry  condition  is  used  at  center  of  the  wire  for 
simulation.  The  outer  surface  of  the  Si  NW  is  assumed  to  be  a 
free  surface,  therefore  the  appropriate  boundary  condition  for 
stress  is 


-"•<r  Rq  =  0  (7) 

At  the  center  of  the  Si  NW,  a  symmetry  boundary  condition  is 
maintained  for  stress. 
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The  outer  surface  of  the  core-shell  structure  is  assumed  to  be  a 
free  surface  without  constraints,  therefore  the  appropriate 
boundary  condition  is 
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At  the  core-shell  interface,  the  radial  stresses  and  displace¬ 
ments  are  taken  to  be  continuous.  At  the  center  of  the  core,  sym¬ 
metry  conditions  are  maintained. 


2.2.4.  Initial  conditions 

The  appropriate  initial  conditions  required  for  simulating  all 
three  case  studies  at  are: 

u\t=Q  =  0  (15) 


2.2.2.  Case  B.  Silicon  nanotube 

For  the  Si  NT  both  the  internal  and  external  radii  defined  as  Ro 
and  Ri  respectively  are  functions  of  time  as  the  structure  expands 
due  to  lithiation.  As  the  inner  surface  of  the  nanotube  is  assumed 
to  be  devoid  of  electrolyte,  any  electrochemical  reactions  at  the 
inner  interface  can  be  ignored  and  therefore  a  zero  flux  condition 
is  used.  At  the  outer  surface,  a  constant  current  flux  condition  is 
used.  Therefore,  the  appropriate  boundary  conditions  are  given 


by: 
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For  the  Si  NT,  both  the  inner  and  outer  surfaces  are  assumed  to 
be  free  surfaces,  therefore  the  proper  boundary  conditions  are  as 


follows: 
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2.2.3.  Case  C.  Si  shell  carbon  core  structures 

For  the  core— shell  structures,  there  are  two  distinct  domains, 
inner  region  representing  the  carbon  core  and  the  outer  region 
representing  the  Si  shell.  The  interface  between  the  core  and 
shell  is  represented  by  R[nt  which  changes  position  with  time.  As 
mentioned  earlier,  we  assume  that  lithium  insertion  inside  the 
carbon  core  is  ignored.  It  should  be  noted  that  in  practical  ap¬ 
plications,  there  is  a  range  of  potentials  in  which  C  and  Si  can  co¬ 
insert  lithium.  However,  for  the  structures  used  in  the  simula¬ 
tions,  the  ratio  of  the  volume  of  carbon  to  Si  is  small  (1/81),  and 
furthermore,  carbon  does  not  undergo  significant  volume 
expansion  (  —  8%)  compared  to  Si.  Therefore,  ignoring  lithium 
insertion  in  carbon  structures  is  a  reasonable  assumption.  The 
model  could  also  be  applicable  to  similar  shell  core  structures  [3]. 
A  constant  current  condition  is  assumed  for  the  outer  surface  of 
the  shell.  Therefore,  the  correct  boundary  conditions  for  outer 
surface  of  Si  shell  and  core/shell  interface  are  as  follows 
respectively: 
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3.  Solving  methodology  &  parameters 

The  model  was  developed  and  executed  in  COMSOL  Multi¬ 
physics.  The  force  balance  equation  (Eq.  (5))  was  solved  in  all 
domains  using  the  in-built  structural  mechanics  module.  The 
material  balance  equation  (Eq.  (2))  was  incorporated  in  a  general 
PDE  interface  to  solve  in  the  respective  domains.  The  1- 
d  problem  was  developed  by  invoking  the  2-d  axisymmetric 
geometry  module  in  COMSOL  Multiphysics  and  imposing  the 
zero  flux  condition  for  all  variables  solved  for  along  the  axial  co¬ 
ordinates.  As  mentioned  in  Part  I  of  the  paper,  the  deformation  of 
the  mesh  is  calculated  directly  from  the  total  displacement 
components  (elastic  and  chemical  strains)  which  are  calculated 
from  structural  mechanics  module,  therefore  making  it  a 
Lagrangian  method.  The  Arbitrary  Lagrangian-Eulerian  method 
was  employed  as  the  mesh  movement  technique.  The  parameters 
used  for  simulation  are  given  in  Table  1.  The  outer  surface  flux 
boundary  condition  has  the  current  density  iapp(t)  which  is 
calculated  based  on  the  ratio  of  total  applied  current  and  the 
varying  outer  surface  area.  The  flux  at  the  outer  surface  is  related 
to  the  current  density  through  Eq.  (6).  For  the  Si  shell  carbon  core 
case,  as  no  lithium  insertion  is  assumed  inside  the  CNT  core,  the 
total  current  is  calculated  based  on  the  initial  volume  of  the  Si 
shell  only.  For  the  Si  NT,  the  initial  volume  of  the  annular  portion 


Table  1 

List  of  parameter  values  used  in  the  simulation. 


Parameter 

Value 

Units 

Diffusion  coefficient  of  Li  in  Si, 
£>us,s  [20] 

2  x  10-12 

cm2  s-1 

Maximum  number  of  moles 
of  Li  that  can  reversibly  alloy 
per  mole  of  Si,  Ax 

15/4 

No  units 

Maximum  insertion  coefficient 
of  Li 

in  LiS,  Azmax 

1 

No  units 

Thermodynamic  factor,  aLis 

1 

No  units 

Partial  molar  volume  of  Li 
in  Si,  VLis 

13.1 

ml  mol-1 

Young’s  modulus,  E  [21,13] 

92.16  (for  LiS,  200 
(for  SWCNT), 

30  (for  MWCNT) 

GPa 

Poisson’s  ratio,  v 

0.27  (for  LiS),  0.35 
(for  SWCNT  &  MWCNT) 

No  units 
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was  considered  for  total  current  calculation.  Single  walled  and 
multi  walled  carbon  nanotubes  were  considered  as  core  mate¬ 
rials  for  Si  shell  carbon  core  structures  which  were  simulated 
with  appropriate  parameters.  These  cases  henceforth  will  be 
addressed  as  SWCNT/Si  and  MWCNT/Si  respectively. 

4.  Results  and  discussion 

In  this  section,  simulation  results  from  the  2-d  model  (presented 
in  part  I)  and  1-d  model  (presented  here)  are  compared.  Subse¬ 
quently,  simulation  results  for  other  nanostructures  are  discussed. 
The  stress  evolution  in  different  type  of  nano-structures  is  also 
discussed. 

4.1.  Electrochemical  lithiation  of  silicon  nanowire  (Si  NW)-l-d  vs  2- 
d  model 

As  stated  earlier,  the  1-d  axisymmetric  plane  strain  model  was 
developed  by  invoking  a  2-d  axisymmetric  geometry  in  COMSOL 
and  imposing  zero  flux  conditions  along  the  axial  co-ordinate 
which  renders  the  1-d  model  valid  only  far  away  from  the  top  and 
base  of  the  2D  geometry.  To  compare  the  results  from  the  1- 
d  model  to  the  2-d  model  (developed  in  part  I  of  the  paper), 
variables  were  evaluated  at  a  radial  cut  at  half  the  height  of  the 
chosen  2-d  axisymmetric  geometry  i.e.,  for  the  2-d  model  the 
stress  profiles  were  obtained  across  radius  at  height  of 
z  =  Hc u  +  Hnw/2  and  were  compared  with  the  results  from  the  1  - 
d  model.  Fig.  2  presents  the  comparison  of  radial  and  tangential 
stress  profiles  obtained  from  1-d  and  2-d  models  at  times  t=  1, 10, 
100  and  1000  s.  The  1-d  model  was  used  to  simulate  the  elec¬ 
trochemical  lithiation  of  Si  NW  of  radius  50  nm  and  height  10p  at 
1C  rate.  From  the  figure,  it  is  evident  that  the  stress  profiles  from 


both  the  1-d  and  2-d  models  agree  well  qualitatively.  Both  models 
predict  the  increase  in  stress  at  start  of  lithiation  and  the  subse¬ 
quent  decrease  over  a  broad  time  scale  on  continuous  lithiation. 
The  tensile  nature  of  the  tangential  and  radial  stress  components 
towards  the  center  and  equality  of  radial  and  tangential  stress  at 
the  center  of  the  wire  are  accurately  captured  by  both  the  models. 
The  compressive  nature  of  the  tangential  stress  at  the  surface  of 
the  nanowire  is  also  captured  by  both  the  models.  At  longer  times, 
the  value  of  radial  and  tangential  stresses  obtained  from  the  1- 
d  model  are  observed  to  be  higher  in  magnitude  than  the  2D 
model.  (-25%  difference  in  radial  and  tangential  stresses  at  the 
center  and  -25%  difference  in  tangential  stresses  at  the  outer 
surface).  This  is  possibly  due  to  the  negligence  of  the  current 
collector  base  in  the  1-d  model.  Also,  the  1-d  plane  strain  model 
ignores  the  effect  of  axial  strain,  which  can  impact  the  radial  and 
the  tangential  components.  It  is  worth  noting  that  the  axial 
stresses  dominate  at  longer  times  (see  part  I),  which  can  be  the 
reason  for  the  increased  deviation  observed  between  1-d  &  2- 
d  models  at  longer  times.  The  deviation  in  stress  profiles  also 
lead  to  differences  in  prediction  of  expanded  radius  (at  the  end  of 
lithiation).  The  increased  radial  stresses  observed  for  the  1- 
d  model  at  longer  times,  results  in  higher  predictions  of  the 
expanded  radius.  There  is  —3.7%  difference  in  prediction  of  the 
expanded  radius  between  the  1-d  and  2-d  models. 

Fig.  3  shows  the  maximum  radial  and  tangential  stress  profiles 
predicted  by  the  1-d  and  2-d  models.  The  structure  simulated  is  Si 
NW  of  radius  50  nm  and  height  10p  at  1C  rate.  The  maximum  radial 
and  tangential  stresses  always  occur  at  the  center  and  outer  surface 
of  Si  NW  respectively.  The  qualitative  features  of  the  maximum 
stress  profiles  e.g.  the  increase  and  decrease  of  stresses  with  time 
and  occurrence  of  peak  stresses  are  in  good  agreement  between  the 
1-d  and  2-d  models.  The  maximum  radial  &  tangential  stresses  are 


u 

tri 


Fig.  2.  Comparison  of  radial  and  tangential  stresses  for  1-d  and  2-d  models  for  Si-NW  at  times:  a)  1  s  b)  10  s  c)  100  s  &  d)  1000  s. 
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(a)  Comparison  of  maximum  radial  stress  profiles  at  cen¬ 
ter  of  Si  NW  (r  =  0)  with  time  for  1-d  and  2-d  models 


Log(t)  (s) 


(b)  Comparison  of  maximum  tangential  stress  profiles  at 
outer  surface  of  Si  NW  (r  =  Rq)  with  time  for  1-d  and 
2-d  models 


Fig.  3.  Plots  comparing  evolution  of  maximum  radial  and  tangential  stress  profiles  with  time  for  Si  NW  from  1-d  and  2-d  models. 


also  in  good  quantitative  agreement,  except  at  longer  times,  i.e., 
towards  end  of  lithiation.  The  possible  reason  for  this  difference 
was  explained  earlier. 

Fig.  4a  and  b  compares  the  maximum  radial  and  tangential 
stresses  obtained  from  the  1-d  and  2-d  models  as  a  function  of  the 
rate  of  lithiation  (C  rate)  of  Si  NW  of  radius  50  nm  and  height  lOp. 
The  C  rates  chosen  were  C/10,  C/2,  C,  2C  and  5C.  Most  of  the  general 
features  such  as  the  increase  and  decrease  of  stresses  with  time, 
increase  in  peak  stresses  and  shift  of  peak  stresses  towards  shorter 
times  with  increasing  rates  are  qualitatively  in  good  agreement 
between  both  the  1-d  and  2-d  models.  The  difference  in  the  radial 
and  tangential  stress  predictions  from  the  models  increases  with 
increase  in  charging  rate  towards  the  end  of  lithiation  (radial  stress 
~1%  difference  &  tangential  stress  ~  1%  difference  at  1C  rate  vs. 
radial  stress  ~  11%  difference  at  5C  rate,  and  tangential  stress  ~  11% 
difference  at  5C  rate  at  t  =  100  s).  This  observation  can  be  explained 
as  an  effect  of  negligence  of  base  and  axial  strain  in  the  1-d  model 
as  mentioned  previously.  Furthermore,  the  tangential  stresses 
predicted  by  the  ID  model  is  higher  than  the  2-d  model,  at  short 


times,  especially  at  higher  rates  (68%  at  1C  rate  vs.  70%  at  5C  rate  at 
t  =  2  ms).  Overall,  as  the  predictions  from  the  1-d  model  compare 
quite  well  to  those  from  the  2-d  model,  therefore,  it  was  used  to 
study  other  nanostructures. 

4.2.  Silicon  nanotube  (Si  NT) 

To  simulate  nanotube  structures,  a  Si  NT  structure  with  an  outer 
radius  of  50  nm  and  an  inner  to  outer  radius  i.e.,  Rf/Ro  of  0.1  was 
considered.  The  simulation  results  from  the  Si  NT  structure  were 
compared  to  the  results  from  Si  NW  (50  nm  radius)  for  lithiation  at 
1C  rate.  Fig.  5  compares  the  radial  stress  profiles  for  Si  NW  to  Si  NT 
at  times  1, 10  and  100  s  respectively.  From  the  plots  it  is  evident  that 
for  the  Si  NT,  the  stresses  are  tensile  along  the  radius  and  tend 
towards  zero  at  the  inner  and  outer  edges  as  a  result  of  free  sur¬ 
faces.  Similar  to  Si  NW,  the  stresses  in  the  Si  NT  increase  at  short 
times  with  initial  lithiation  (<10  s),  and  decrease  at  longer  times 
(>10  s)  towards  the  end  of  lithiation.  For  the  Si  NW,  the  maximum 
radial  stress  always  occurs  at  the  center,  whereas  for  the  Si  NT,  the 
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(a)  Comparison  of  maximum  radial  stress  profiles  at  cen-  (b)  Comparison  of  maximum  tangential  stress  profiles  at 
ter  of  Si  NW  (r  =  0)  with  time  for  1-d  and  2-d  models  outer  radius  of  Si  NW  (r  =  R0)  with  time  for  1-d  and 

as  a  function  of  different  rates  of  lithiation  2-d  models  as  a  function  of  different  rates  of  lithiation 


Fig.  4.  Plots  comparing  evolution  of  maximum  radial  and  tangential  stress  profiles  of  Si  NW  with  time  for  1-d  and  2-d  models  as  a  function  of  different  rates  of  lithiation. 
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Radial  Position  (nm) 


Fig.  5.  Comparison  of  radial  stress  profiles  for  Si  NW  and  Si  NT  at  various  times  t=  1,10  Hg-  7.  Comparison  of  tangential  stress  profiles  for  Si  NW  and  Si  NT  at  various  times 

and  100  s.  t=l,  10  and  100  s. 


maximum  radial  stress  occurs  somewhere  along  the  radius  as  a 
function  of  time  and  moves  closer  to  the  inner  radius  with 
increasing  time.  Also,  the  local  radial  stresses  developed  across  the 
radius  in  the  Si  NT  is  in  general  lower  than  that  of  the  stresses  in  the 
NW.  A  comparison  of  maximum  radial  stress  profiles  of  Si  NW  and 
Si  NT  with  time  is  presented  in  Fig.  6.  Furthermore,  Fig.  5  also  shows 
that  lithiation  of  the  Si  NT  causes  the  inner  interface  to  move  away 
from  the  initial  position  (not  obvious  from  the  plot  in  Fig.  5).  At 
t  =  100  s,  the  inner  surface  has  a  positive  displacement  of  about 
0.102  nm  from  its  initial  position. 

Fig.  7  compares  the  tangential  stress  profiles  across  the  radius  of 
the  Si  NW  and  Si  NT  at  times  1, 10  and  100  s.  The  tangential  stresses 
in  Si  NT  are  tensile  at  the  inner  surface  and  compressive  at  the  outer 
surface,  similar  to  the  behavior  of  Si  NW.  The  tangential  stress  in  Si 
NT  structure  is  higher  in  magnitude  closer  to  the  inner  radius 
compared  to  the  tangential  stress  at  center  of  the  Si  NW.  In  fact,  at 
the  interface,  the  Si  NT  tangential  stress  is  significantly  higher  than 


the  Si  NW  center  tangential  stress,  however  closer  to  the  surface, 
the  Si  NT  tangential  stress  is  lower  in  magnitude  compared  to  the  Si 
NW. 

The  1-d  model  was  used  to  study  the  effect  of  Ri/Ro  i.e.,  inner  to 
outer  radius  ratio  on  the  radial  stress  profiles  for  the  Si  NT  of  fixed 
initial  outer  radius  of  50  nm  for  1C  rate.  The  radius  ratios  consid¬ 
ered  were  0.05,  0.1,  0.2,  0.3  and  0.4,  which  correspond  to  volume 
ratios  (void  volume  to  Si  volume)  of  0.0025,  0.01,  0.042,  0.099  and 
0.19  respectively.  The  radial  stresses  at  t  =  10  s  are  shown  in  Fig.  8. 
For  each  case,  the  local  radial  stresses  developed  across  the  radius 
in  the  Si  NT  are  lower  in  magnitude  for  larger  radius  ratios. 
Furthermore,  the  peak  radial  stress  shifts  away  from  the  inner 
surface  for  Si  NT  structures  with  larger  radius  ratios. 

Based  on  the  simulation  results  of  the  Si  NT,  the  radial  stresses 
are  lower  compared  to  the  Si  NW,  while  the  tangential  stresses 
(interface  and  regions  close  to  interface)  are  higher  in  Si  NT 
compared  to  Si  NW.  Also,  the  radial  stresses  are  lower  in  magnitude 


Fig.  8.  Comparison  of  radial  stress  profiles  for  Si  NTs  of  inner  to  outer  radius  ratio  /?,•/ 
Fig.  6.  Comparison  of  maximum  radial  stress  profiles  for  Si  NW  and  Si  NT  with  time.  R0  =  0.05,  0.1,  0.2,  0.3  and  0.4  at  time  t  =  10  s. 
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for  increasing  Ri/Ro  values  for  Si  NT.  Consequently,  more  simulation 
runs  were  performed  to  compare  radial  and  tangential  stresses  of  Si 
NT  to  Si  NW  for  various  Ri/Ro  values.  Fig.  9  shows  the  magnitudes  of 
maximum  tensile  stresses  for  the  Si  NT  for  different  Ri/Ro  ratios.  The 
initial  outer  radius  for  both  NW  and  NT  was  50  nm.  For  the  Si  NT, 
the  tangential  component  of  the  maximum  tensile  stress  occurs  at 
the  inner  surface  while  the  maximum  radial  stress  component 
occurs  somewhere  along  the  radius.  The  maximum  tensile  stress 
for  the  Si  NW  always  occurs  at  the  center  where  the  radial  and 
tangential  stresses  are  equal  in  magnitude.  As  mentioned  earlier, 
the  radial  component  of  the  tensile  stress  decreases  with  increasing 
inner  to  outer  radius  ratios  for  the  Si  NT.  Interestingly,  the 
tangential  component  of  the  maximum  tensile  stress  also  follows  a 
similar  trend.  For  the  case  when  the  radius  ratio  becomes  zero  i.e., 
Ri/Ro  =  0  which  represents  Si  NW,  the  radial  and  tangential  com¬ 
ponents  of  the  maximum  tensile  stress  become  equal  as  a  result  of 
the  symmetry  boundary  condition  represented  by  the  solid  hori¬ 
zontal  line  .In  general,  although  the  radial  component  of  the 
maximum  tensile  stress  are  always  smaller  in  the  case  of  Si  NT 
compared  to  Si  NW,  the  tangential  component  is  larger  than  that  of 
Si  NW.  This  scenario  reverses  for  radius  ratio  values  greater  than 
0.4,  where  both  the  radial  and  tangential  components  of  maximum 
tensile  stresses  for  Si  NT  are  smaller  than  Si  NW.  Such  a  design  is 
beneficial  to  minimize  the  overall  stress  induced  in  nano-structure. 
However,  increase  in  the  hollow  space  for  the  Si  NT  of  fixed  initial 
outer  radius  decreases  the  charge  capacity.  For  a  radius  ratio  of  0.4, 
the  volumetric  capacity  of  Si  NT  decreases  by  16%  compared  to  Si 
NW. 

4.3.  Silicon  shell  carbon  core  structure 

The  1-d  model  was  also  used  to  study  core-shell  structures.  Use 
of  core-shell  carbon/Si  nanostructures  for  anodes  for  use  as 
negative  electrodes  in  lithium-ion  batteries  have  been  reported  in 
literature  [11,12,9].  Cui  and  co-workers  [12]  reported  synthesis  of 
amorphous  silicon  using  chemical  vapor  deposition  technique  onto 
carbon  nano-fibers,  casted  on  stainless  steel  substrates  to  generate 
carbon/Si  core-shell  structures.  Kumta  et  al.  [9]  reported  synthesis 
of  vertically  aligned  Si/carbon  nanotube  arrays  directly  on  Inconel 


Fig.  9.  Comparison  of  maximum  tensile  stresses  of  Si  NT  with  increasing  K//i?0  ratios. 
The  thick  red  line  represents  the  maximum  tensile  stress  at  the  center  of  Si  NW  (/?,■/ 
R0  =  0).  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is 
referred  to  the  web  version  of  this  article.) 


metal  current  collector,  using  a  two-step  chemical  vapor  deposition 
process.  Justin  et  al.  [11]  reported  synthesis  of  smooth  and 
conformal  thin  films  on  various  multi-walled  carbon  nanotube 
using  radio  frequency  plasma  enhanced  chemical  vapor  deposition. 
For  the  simulations  considered  here,  single-walled  carbon  nano¬ 
tube  (SWCNT)  and  multi- walled  carbon  nanotube  (MWCNT) 
structures  were  chosen  as  the  materials  for  core,  and  amorphous  Si 
as  the  shell  material.  The  Young’s  modulus  of  SWCNT  and  MWCNT 
are  200  and  30  GPa  respectively  [13].  The  SWCNT(MWCNT)/Si 
core-shell  structure  has  an  outer  radius  of  50  nm  with  a  SWCNT 
(MWCNT)  core  of  radius  R[nt  =  1/10%  be.,  5  nm.  Fig.  10  plots  the 
radial  stress  profiles  for  SWCNT/Si  and  MWCNT/Si  structures  at 
times  of  1, 10,  and  100  s,  respectively.  For  both  these  structures,  the 
radial  stresses  are  always  maximum  at  the  center  extending  up  to 
CNT/Si  interface.  This  is  similar  to  Si  NW  structure,  where  the 
maximum  radial  stress  always  occurs  at  the  center,  however  this 
maximum  value  extends  up  to  the  CNT/Si  interface  because  of  the 
constant  stresses  within  the  core  (Appendix  I).  The  maximum 
observed  radial  stress  at  the  CNT/Si  interface  extends  into  the  core, 
and  does  not  vary  spatially  within  the  core.  Since  the  core  does  not 
insert  lithium,  therefore  under  such  circumstances,  the  radial 
stresses  should  remain  constant  based  on  the  force  balance  equa¬ 
tions  (see  Appendix  I).  Furthermore,  the  radial  stresses  in  the  core 
and  the  shell,  increases  sharply  with  time  which  is  quite  different  in 
behavior  from  other  nanostructures,  Si  NW  and  Si  NT.  This  is  due  to 
the  discontinuity  of  concentration  variable  at  the  interface  com¬ 
bined  with  the  absence  of  free  surface.  The  magnitude  of  radial 
stresses  are  larger  for  SWCNT/Si  structure  compared  to  the 
MWCNT/Si  structure  due  to  the  lower  magnitude  of  Young’s 
modulus  of  the  MWCNT  compared  to  SWCNT.  The  absence  of 
chemical  strain  and  the  high  value  of  elastic  moduli,  resists 
movement  of  the  interface.  The  interface  does  not  move  up  to  first 
100  s  of  lithiation. 

Fig.  11  shows  maximum  radial  stress  evolution  of 
SWCNT(MWCNT)/Si  structures  with  time  and  compares  them  to  Si 
NW.  The  maximum  radial  stress  magnitudes  for  both  Si  shell  car¬ 
bon  core  structures  are  much  higher  compared  to  Si  NW.  It  is  to  be 
noted  that  at  very  short  times  (<2  s),  MWCNT/Si  maximum  radial 
stress  is  lower  than  Si  NW  maximum  stress.  This  is  because  at  those 
short  times,  chemical  strain  component  is  almost  insignificant.  So 
the  strains  are  mainly  elastic  strains.  Since  MWCNT  has  the  lowest 
elastic  modulus  compared  to  other  materials  considered  here, 


Fig.  10.  Comparison  of  radial  stress  profiles  for  Si  SWCNT  and  Si  MWCNT  at  various 
times  t  =  1, 10  and  100  s. 
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Fig.  11.  Comparison  of  maximum  radial  stress  profile  for  Si  NW  and  core  radial  stress 
profiles  for  Si  SWCNT  and  Si  MWCNT  with  time. 

therefore  it  generates  lower  stress.  Fig.  12  compares  the  tangential 
stress  profiles  for  SWCNT/Si  and  MWCNT/Si  structures  at  times  of  1 
and  10  s.  For  both  the  core  shell  structures,  tangential  stresses  are 
tensile  inside  the  core  which  increase  with  time  while  at  outer 
surface  they  are  fully  compressive.  This  behavior  is  similar  to  Si  NW 
structure.  Much  higher  stresses  are  observed  as  a  result  of  the  non¬ 
insertion  core.  The  tangential  stress  for  both  structures  become 
discontinuous  at  the  core— shell  interface  and  turns  more 
compressive  with  time.  The  tangential  stresses  for  the  MWCNT/Si 
structure  are  less  compared  to  the  SWCNT/Si  structure.  For 
example,  at  t  =  100  s,  the  maximum  tangential  stress  for  the 
MWCNT/Si  structure  has  a  magnitude  of  612  MPa  compared  to 
1024  MPa  for  the  SWCNT/Si  structure. 

Fig.  13  shows  a  summary  plot  comparing  the  radial  stress  pro¬ 
files  for  different  structures  discussed  earlier.  The  plots  suggest  that 
core  shell  structures  with  non-insertion  core  clearly  develop  higher 
stresses  compared  to  Si  NT  and  Si  NW.  On  the  other  hand,  Si  NT 
structure  generates  lower  stresses  than  Si  NW  and  CNT/Si  core- 


Fig.  12.  Comparison  of  tangential  stress  profiles  for  Si  SWCNT  and  Si  MWCNT  at  times 
t  =  1  and  10  s. 


Radial  Position  (nm) 

Fig.  13.  Comparison  of  radial  stress  profiles  for  Si  NW,  Si  NT,  Si  SWCNT  and  Si  MWCNT 
at  times  t  =  10  s. 

shell  structures  (for  similar  current  densities)  which  could  be  a 
preferred  structure  for  practical  use.  Core-shell  structures  with 
substantial  lithium  insertion  cores  will  behave  differently  from  the 
structures  considered  here  because  of  the  continuous  lithium 
concentration  gradients  across  the  interfaces.  This  will  be  a  subject 
of  future  work. 


5.  Conclusion 

A  1-d  axisymmetric  plane  strain  stress  model  coupled  with 
pressure  induced  diffusion  was  developed  to  analyze  stress  distri¬ 
bution  during  lithium  insertion  in  different  silicon  based  nano¬ 
structures  (Si  NW,  Si  NT,  Si  shell  carbon  core).  The  results  from 
simulation  of  lithiation  of  Si  NW  using  the  1-d  model  are  compared 
to  the  results  from  the  2-d  model  presented  in  part  I  of  the  paper. 
The  two  models  agree  very  well  at  short  times,  however  tend  to 
deviate  slightly  at  longer  times.  Also,  the  modeling  results  infer  that 
Si  NT  could  be  a  better  structure  in  terms  of  mechanical  stability  as 
it  generates  lower  stresses  compared  to  Si  NW  and  CNT/Si  core¬ 
shell  structures.  Furthermore,  core— shell  structures  (SWCNT/Si  and 
MWCNT/Si)  exhibit  very  high  stresses  at  the  core— shell  interface. 
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Appendix  I 

The  radial  and  tangential  stresses  inside  the  core  for  a  Si  shell 
carbon  core  structure,  are  both  constant  across  radius.  Moreover,  in 
the  core  the  radial  and  tangential  stresses  are  equal.  In  this  section, 
a  mathematical  approach  is  explained  to  verify  this  behavior.  In  the 
carbon  core,  Eq.  (5)  is  valid.  In  addition,  two  other  assumptions  are 
made  here: 

•  As  only  elastic  strains  exist  in  non-insertion  carbon  core,  and  the 
deformation  is  very  small,  linear  theory  of  elasticity  is  assumed 
to  hold  good. 

•  Movement  of  core— shell  interface  is  neglected. 
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Replacing  the  stress  terms  in  Eq.  (5)  with  respect  to  displace¬ 
ment  u  we  get 


->d2u 


du 

-f -r-= —  u  =  0 


dr2  dr 


(17) 


The  general  solution  of  this  kind  of  ODE  is  u  =  C\r  +  C2I /r  where 
Ci  and  C2  are  constants  of  integration.  Applying  symmetry  bound¬ 
ary  condition  at  r  =  0,  the  final  solution  is  u  =  C\r.  Plugging  this 
solution  back  to  the  radial  and  tangential  stress  expressions,  it  can 
be  proved  that  orr  and  goo  are  constant  and  equal  to  each  other  in 
the  core. 


Appendix  II 

For  development  of  both  2-d  and  1-d  model,  we  split  the  total 
strain  eT  into  chemical  and  elastic  components  i.e.,  eCh  and  ee\.  This 
methodology  is  extensively  reported  in  literature  [14—18  .  Chris¬ 
tensen  et  al.  [19  reports  a  slightly  different  technique,  where  in  the 
total  stress  is  split  into  an  elastic  component  and  thermodynamic 
pressure.  The  thermodynamic  pressure  p  is  solved  as  a  separate 
variable.  The  author  claims  that  this  approach  is  more  general  as  p 
is  not  restricted  to  being  just  a  function  of  stress  tensor.  In  this 
work,  the  thermodynamic  pressure  p  is  assumed  to  be  equivalent  to 
the  mean  normal  pressure  which  is  function  of  stress  tensor  as 
discussed  earlier  in  Part  I.  To  understand  the  difference  between 
either  of  the  approaches,  we  developed  the  1-d  axisymmetric 
model  based  on  the  procedure  of  splitting  the  stress,  and  compared 
to  the  alternate  approach  of  splitting  the  strain.  Fig.  14  shows  the 
comparisons  of  maximum  radial  and  tangential  stress  profiles  with 
time  of  Si  NW  of  radius  50  nm  from  both  the  1-d  models  developed 
following  the  stress  and  strain  splitting  procedures.  Both  ap¬ 
proaches  show  very  good  qualitative  comparison  -  maximum 
radial  stress  at  center  and  maximum  tangential  stress  variation  at 
outer  surface  is  predicted  by  either  approach.  At  short  times,  the 
radial  stresses  predicted  by  the  stress  splitting  method  is  slightly 
higher  than  the  strain  splitting  method,  while  the  tangential 
stresses  are  quantitatively  very  close.  The  maxima  (with  time) 
predicted  for  both  radial  and  tangential  stresses  using  the  strain 
splitting  approach  is  larger  compared  to  the  stress  splitting 
approach.  Moreover,  the  occurrence  of  the  peak  stresses  is  also 
shifted  to  the  right  for  the  strain  splitting  method.  The  pressure 
gradient  in  the  Li  flux  equation  is  a  dominating  factor  controlling 


Log  (t)  (s) 


Fig.  14.  Comparison  of  maximum  radial  and  tangential  stress  profiles  of  Si  NW  from  1  - 
d  models  (strain  and  stress  splitting). 


lithium  insertion  which  affects  the  thermodynamic  pressure  and  in 
turn  the  stresses.  This  effect  becomes  prominent  when  the  Li 
concentration  gradients  are  high  at  shorter  times  which  explains 
the  deviation  between  the  two  approaches.  At  prolonged  times,  the 
concentration  gradients  become  more  uniform  which  prevents  the 
profiles  from  diverging  at  longer  times. 

Nomenclature 

Ct  total  concentration  of  species,  mol  m~3 

N  flux  of  species  S,  or  LiS,  mol/m2s 

Xus  mole  fraction  of  species,  LiS 
A z  insertion  coefficient  of  Li  in  LiAZSii/A* 

Ax  maximum  number  of  moles  of  Li  that  can  reversibly  alloy 
per  mole  of  Si 
t  time,  s 

r  radial  co-ordinate,  m 

6  tangential  co-ordinate,  deg 

z  axial  co-ordinates,  m 

Dus.s  binary  diffusion  coefficient  of  Li  in  LiS,  m2  s-1 
aus  thermodynamic  factor 

Mys  molar  mass  of  species,  LiS,  g  mol-1 
F  Faraday’s  constant,  C  g-1  equiv 

R  universal  gas  constant,  J  mol-1  I< 

T  temperature,  I< 

Si  stoichiometric  coefficient 

VLis  partial  molar  volume  of  species,  LiS,  m3  mol-1 
p  density  of  LizSii/Ax,  g  rn-3 

p  thermodynamic  pressure,  N  m-2 

u  displacement  vector  field 

u  displacement  in  r  co-ordinate,  m 

n  unit  outward  normal  vector,  m 

iapp  current  density  related  to  the  outer  surface  area  of  the 
nanostructures,  A  m-2 
t  stress  tensor  matrix 

a  stress  components  of  the  stress  tensor  matrix,  N  m-2 

e  strain  tensor  matrix 

e  strain  components  of  the  strain  tensor  matrix,  N  m-2 

X  lame  constant  given  by  Ev\\  +  v,  N  m~2 

p  Lame  constant  given  by  E/2(l  +  *0(1  -  2z0,  N  m“2 

E  Young’s  modulus,  N  m-2 

v  Poisson’s  ratio 

Subscripts 

LiS  lithiated  host 

S  unoccupied  host 

T  total 

max  maximum 
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